Density response of a trapped Fermi gas: a crossover from the pair vibration mode to 

the Goldstone mode 



A. Korolyuk^, J. J. Kinnunen^, P. Torma^-^ 
^Department of Applied Physics, School of Science, 
Aalto University, P.O. Box 15100, 00076 Aalto, Finland 
^ Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106-4030, USA 

We consider the density response of a trapped two-component Fermi gas. Combining the 
Bogoliubov-deGennes method with the random phase approximation allows the study of both col- 
lective and single particle excitations. Calculating the density response across a wide range of 
interactions, we observe a crossover from a weakly interacting pair vibration mode to a strongly 
interacting Goldstone mode. The crossover is associated with a depressed collective mode fre- 
quency and an increased damping rate, in agreement with density response experiments performed 
in strongly interacting atomic gases. 
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I. INTRODUCTION 

The response of a many-body system to external per- 
turbations can be understood in terms of collective, i.e. 
many-body, and single particle excitations. Collective 
modes of a quantum fluid are well described by the hy- 
drodynamic model in the strongly interacting hydrody- 
namic and the weakly interacting collisionless limits [ij- 
[^. However, experimental studies on collective modes of 
trapped Fermi gases yielded surprises [l0l - [l2| that 
did not fit in the simple picture obtained from the hydro- 
dynamic theory. The difference was suggested to lie in 
the interplay between the single particle excitations and 
the collective modes [Tsj . 

The Bogoliubov-deGennes (BdG) mean-fleld theory 
provides a microscopic description of the trapped Fermi 
gas. While the single particle excitations are readily ac- 
cessible from the basic BdG theory, also the collective 
modes can be obtained by using the random phase ap- 
proximation [l^ (RPA). The method allows the study of 
different modes p^l - [l8| . but here we will concentrate on 
the monopole mode in a spherically symmetric trapped 
Fermi gas. The method has already been used for study- 
ing the monopole mode of a trapped Fermi gas in the lim- 
its of weak and strong interactions. The purpose of the 
present work is to study the interesting crossover region 
between the weakly interacting and the strongly interact- 
ing regimes where the hydrodynamic theory has failed to 
properly describe the experiments. 

The weakly interacting limit was studied in Ref. [l9| . 
and in this limit the lowest lying collective mode with the 
monopole symmetry was identified as the pair vibration 
mode with the frequency u! = 2A{0)/H, where A(0) is 
the superfluid excitation gap in the center of the trap. 
This mode was argued to be the precursor of the Gold- 
stone mode in the strongly interacting regime with the 
frequency 2 lot , where wt is the harmonic trapping fre- 
quency, in agreement with the strongly interacting hy- 



drodynamic theory. The Goldstone mode in this more 
strongly interacting regime was considered in Ref. [20j . 

Here we study the actual transition between the two 
regimes and how the collective pair vibration mode trans- 
forms into the collective Goldstone mode when the in- 
teraction strength is increased. The crossover region 
is characterized by a depression of the collective mode 
frequencies and increased damping, in qualitative agree- 
ment with experiments done in non-spheric ally sym met- 
ric traps and with various collective modes |10l - [l3 |. The 
full crossover could be studied in future experiments, and 
the study of the lately realized systems with small num- 
bers of atoms [2l| might be helpful in finding the required 
parameters. 

This paper is organised as follows. In Section [H] we 
review the standard theory of linear response for a small 
disturbance to the system and discuss the connection be- 
tween the density response and the frequency of collective 
excitations. In Section IIIII we discuss the details of the 
methods used. In Section ITVl we present and analyse our 
results and Section |V] summarises the main findings and 
conclusions of this work. 



II. LINEAR RESPONSE AND COLLECTIVE 
BEHAVIOUR 

In this section we review the calculation showing the 
connection between the density response and the fre- 
quency of collective excitations of a system. This shows 
how peaks in the imaginary part of the density response 
function correspond to the frequencies of collective exci- 
tations, providing the theoretical background for inter- 
preting the results of the numerically calculated density 
response function. 



A. The density response function 
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Let US consider a system described by a full many- 
body Hamiltonian H. It is assumed that the system is 
initially in a pure state \4>o) which is an eigenstate of 
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the Hamiltonian H , for example the ground state. At 
the moment t — Q a, probing field or a pertm'bation V 
is switched on. From that moment the system evolves 
under the Hamiltonian H' = H + V . We denote the state 
of the system at the time t as \(t>{t)) ■ In the absence of 
any probing field, the state of the system at time t would 
remain |(/)o). The difference between \<t>{t)) and |(/)o) is 
thus caused by the probing field V . This difference can be 
measured with the help of an observable 0(r). Without 
the probing field, the average of the operator 0(r) (the 
statistical average after a series of measurements) would 
be 



Oo(r) = (<^o|O(r)|0o) 



(1) 



When the field V is switched on, the measurements would 
yield 



The difference 



50{v,t)^0{v,t)~Oo{v) 



(2) 



(3) 



indicates how the perturbation V has influenced the 
physical observable O. 

The difference, or response, 50 can be calculated in 
the interaction picture representation. In the density re- 
sponse approach we consider as a weak disturbance to 
the system. In this case we can take into account only 
the first (linear) order of the perturbation V . Thus the 
response is 



50{v,t)^-i / dt'iU 6i{v,t),Vi{t') I0o), (4) 



where Oi and Vi are the operators O and V in the inter- 
action picture representation. 

In this work we are interested in the particular case of 
an external potential w(r,t) that couples to the density 
/5(r), corresponding to an operator V = J drp{r)v{r,t). 
Then SO becomes 
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SO{r,t)= J dt'dr'Air,r',t,t')v{r',t'), 



(5) 



where the kernel 



A{r,r',t,t')=-i{(bo\ [Oi(r,i),Pi(r',t')J \MHt-t') 

(6) 

is called the response function. 

The response function A{r, r' ,t, t') in Eq. ^ shows the 
change in the observable O measured at the point (r, t) 
due to the infinitesimally small perturbation at the point 
(r',i'). If the Hamiltonian H does not depend on time, 
then the time-dependence of A will be on the difference 
t - t' only. 



The derivation above holds for any general observable 
O. In the special case O — p, i.e. when the observed 
quantity is also the density, the response function is called 
the density response function, and the expression for it 
is 



A{r,r',t~t') 



-i{4,o\[pi{r,t),pi{r',t')]\c^o)0{t~t') 



(7) 



B. The frequency of the collective excitations in 
the density response 

After a Fourier transformation of Eq. ([7]) one obtains 



A{r,r\uj)^~i j die-'"*(0o|[pi(r,O),pi(r',t)]|0o). 

— OO 

We assume, formally, that the Hamiltonian H is diag- 
onalized, and En and \n) are its eigenvalues and eigen- 
vectors: H\n) = En\n), n = 0, 1, As eigenstates 

of a hermitian operator, the eigenvectors are orthogonal 
(n| n') = 5n,n' and form a complete basis \n) (n| = 1. 
One can assume that the initial state of the system, ear- 
lier denoted as |0o), is the ground state |0): \(f)o) = |0). 

Transforming back to the Schrodinger picture repre- 
sentation pi{r,t) — e*^*/5e~*^*, Eq. ([SJ becomes 



^(r,r',a;) = 2 



dt e 



(9) 



Im(e-'^°*(</>o|p(r)e^^*/}(r')|0o)) 



Using the completeness of the eigenbasis and performing 
the time integral yields 



A{r,r',u}) = 



^ - {En - Eo) 

{Mp{r')\n) (n|p(r)|0o) 
uj + {En- Eo) 



(10) 



Thus, the density response ^(r, r',w) yields informa- 
tion about the excitation spectrum of the system En—E^. 
In Section IIVI we will analyse the numerically calcu- 
lated density response function and identify the poles of 
Eq. (fTO|). i.e. the frequencies w in which A has peaks, 
as the collective or single particle excitation frequencies. 
Note that here, formally. En includes all excitations of 
the system, both collective and single particle ones. 



III. THE RANDOM PHASE APPROXIMATION 
FOR THE HARMONIC TRAP GEOMETRY 

In the previous section we discussed the general defini- 
tions of density response function and its connection to 
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the collective excitation frequencies. Here we will apply 
the random phase approximation (RPA) HI] to ex- 
press the density response using Green's functions. These 
in turn will be solved using the Bogoliubov-deGennes 
(BdG) mean-field theory, allowing an efficient solution of 
the density response for a spherically symmetric trapped 
atomic gas. 



Density response expressed using Green's 
functions 



A dilute ultracold two-component atomic gas trapped 
in a spherically symmetric trap can be described with the 
Hamiltonian 



H = K + U, 



where 
K = 



E 



drV'i(r) 



2m 



(11) 



Mr) (12) 



is the single particle Hamiltonian containing the kinetic 
energy and the harmonic trapping potential of frequency 
'ipa\T^) is the annihilation (creation) operator of an 
atom in a state a at point r, and 



U = 



\ j rfrdrVl(r')^;(r')V'Mr)„V(r)5 (r' - r) 



(13) 

describes two-body interactions. The short-range in- 
teractions can be described by the two-body scatter- 
ing T-matrix in the contact potential approximation 
g (r' — r) = g^S (r' — r), where the coupling constant g^ 
is related to the physical s-wave scattering length through 
the relation 1/go = ^ - Efc 2^ l2l- Eq. dH]) becomes 



E dvi^i{v)4{v)Mr)Mr)- (14) 



For the calculation of the density response we need to 
add a weak probing field V (for more details see IH Ap . 
We consider 



V 



dr [4)^{r, t)n^{r) + (/«;(r, t)n^(r) 
+?7(r,i)Vi(r)Vt(r) + ff.c.]. 



(15) 



The three terms in V allow the calculation of 
the density responses to the three external fields 
(fi^ir, t), 04,(r, t),r]{r, t), which couple to the system in dif- 
ferent ways. To simplify the notation, we consider a re- 
sponse to a generic field h{r,t), where h{r,t) can be any 
of these three external fields. 

As physical observables for which we will analyse the 
response, we consider the densities and pi of up and 



down components, respectively. Additionally, we con- 
sider the response of the order parameter A(r) to the 
probing field. We expect that the same frequencies of 
collective excitations will appear in all three responses; 
naturally, this can be confirmed numerically. 

A useful way of formulating the density response is 
by using the two-point Green's functions Gij(l,2) = 

— (T\I'i(l)4'J(2)\ in the Nambu formalism, where 1 de- 



notes space-time point xi,ti and ^'(1) = 



^t(l) 
V'I(l) 



. The 



Green's function for 1 = 2 has a simple physical meaning: 



G(l,l) = 



_f Ptil) A(l) 



A*(l) -Piil) 



(16) 



where p-^{l), p^{l) are the densities and A(l) is the order 
parameter, i.e. the gap at the point 1 = (xi,ti). 
We consider the response of the Green's function: 



-t-00 



5G,y (l,2)= J dhj dx5A,,(l,2,5)/i(5). (17) 



This can be interpreted as the change of the Green's func- 
tion, or how the transition between the points (xi,ti) 
and (x2,t2) alters due to a probing field h at the point 
(x5,t5). The kernel can be written as 



A,, (1,2, 5) 



5Gy(l,2) 



Sh{5) 



(18) 



For the one-point Green's functions the response is 



+00 

5G,,(1,1)= J dt5 J^dx5Aj{l,5)h{5), 



where 



and 



A,,(1,3) = A,,(1,1,3) 
<5Gy-(l,l) 



A,, (1,5) 



6h{5) 



(19) 



(20) 



(21) 



As in Eq. (|16l) . this is a 2x2 matrix with the differ- 
ent elements describing the response to different fields. 
For example, the response of the density of up particles 
A^^{1, 5) is 



+00 



%(1)= / dU / d^5Anil.5)hi5). (22) 



Starting from the Hamiltonian (ITT|) and using the 
random-phase approximation (RPA) |22|, we write (for 
more details see Appendix 
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^,(1,2,5) = Ao.j(l,2,5)+go5I / d3L,fcfcj (1, 2, 3)Ai,(3, 3, 5) - 50 / d3L,fci,(l, 2, 3)^^,(3, 3, 5), (23) 

fc,; •' k,i •' 

where L,fc;j(l,2,3) = G,fc(l, 3)G;j(3, 2), and Ao(l, 2, 5)^= L^ttj (1> 2, 5), £,^^^(1, 2, 5), i,tij(l, 2, 5) + 2, 5) 

for /i = 0^, ?7 respectively. For the density response Aij{l, 5) the equation in Eq. (j23p becomes 




where Life/j(l,3) = Gife(l, 3)G;j(3, 1) and Aojj(1,5) is 
correspondingly expressed via Likij{l, 3). 

B. Spherical symmetry and the angular momentum 

Solving even the mean-field density response in 
Eq. (IM)) in a spatially inhomogeneous system such as the 
trapped gas is a numerically very demanding task. How- 
ever, taking advantage of the assumed spherical symme- 
try of the system, the Eq. ([24]) can be greatly simplified 
by expressing all quantities in spherical coordinates. The 
spherical symmetry implies separation of different angu- 
lar momentum responses, allowing a significant simplifi- 
cation of the numerical calculations. 

Expressing the response function using Legendre poly- 
nomials (and using the knowledge that the response func- 
tion depends only on the time difference ti — 12) yields 

i,fe (1, 2) = ^ Afc,L (r-i , r2 , uj)Pl (cos ^)e-''^^'' -''\ 

L.UJ 

(25) 

where w is a frequency (can be either real or imaginary) , 
Pl (cos 7) is the Legendre polynomial, and 7 is the angle 
between the vectors ri and r2. Other functions involved 
in the solution, such as the Green's function Gifc(l,3) 
and iifcij(l, 2, 3) can be decomposed in the same way. 
Instead of eight parameters ri,ti,r2,t2 all the functions 
now depend only on four parameters: ri,r2 (magnitudes 
of the vectors ri and r2), w (frequency) and L (angular 
momentum). 

With this decomposition, Eq. ([25)1 simplifies to 

+ go ^ j rldr3Cikkjx{ri,r3,uj)Aii,L{r3,r5,uj) 

~ 90 ^ J rldrsdkij^L {ri ,r3,uj)Aki,L{r3,r5,uj). 

(26) 

As explained in Section fllBl the peaks in the density 
response ^ij,L(?'i,?'5,w) yield the frequencies of collec- 
tive excitations, now corresponding to a specific angular 
momentum L. For example, in this work we will study 
the case L = 0, in other words the spherically symmetric 
monopole mode collective excitation. 



In the Section IIII CI we calculate the coefficient 
'Cifcfcj,L(ri, r3, w) via the Bogoliubov-deGennes approxi- 
mation. In Section IIIIDI we consider the simplifications 
we need to proceed from the analytical equations to the 
numerical calculations. 



C. Green's functions in the Bogoliubov-deGennes 
approach 

In this section we will use the mean-field Bogoliubov- 
deGennes (BdG) method to obtain the single particle 
Green's functions in the trap, thus allowing the numer- 
ical solution of the density response. The Bogoliubov- 
deGennes equations can be obtained by approximat- 
ing the interaction part of the Hamiltonian (I13|) by a 
quadratic form 

U = - J dr V'j(r)V'|(r)A (r) + H.c, (27) 

where A (r) = \go\ {'i(j^{r)^i{r)) is the or- 
der parameter. The Hartree energy is 

5o/rfr (^?/'j(r)V't(r)'^i(r) +V'|(r)-0i(r)'^t(r))> where 
n-|-(_l,)(r) is the density of up(down) particles. Most of 
the results below neglect the Hartree energy because, 
especially for strong interactions \k-pa\ ~ 1, the inclusion 
of the Hartree energy would require a prohibitively 
high cutoff energy, which is a critical quantity for the 
efficiency of the numerical method. However, we will 
briefly consider below the effect of the Hartree energy 
shift in the weakly interacting regime. 

The present work is concentrated in the moderate to 
weak interaction regime \k-pa\ < 1, and the effect of the 
Hartree energy is to compress the gas slightly, increas- 
ing the order parameter at the centre of the trap [23 |. 
As will be seen in Section IIVI the key quantity for the 
density response is the order parameter proflle instead of 
the interaction strength. Consequently the results will be 
shown as a function of the order parameter at the cen- 
tre of the trap A(r — 0). The compressing effect of the 
Hartree energy is thus relevant for finding the precise cor- 
respondence between the system input parameters (atom 
numbers, masses, trap frequencies, temperature, and in- 
teraction strength) and the order parameter. However, 
as will be seen, the Hartree energy does not change the 
qualitative features of the density response. 
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With these approximations the Hamiltonian in 
Eq. pT|) becomes quadratic in V'a(r): 



2m 



dr V'|(r)V'|(r)A(r) +i?.c. 



(28) 



The fermion fields can be expressed in the harmonic os- 
cillator basis as = J2nim Rni{r)Yi^{9)cnima, where 
the operator Cnima destroys an atom from the harmonic 
oscillator eigenstate nlm, Yim{0) are the spherical har- 
monics, and the radial eigenstates are given by 



Rni{r)^V2{mL0Tf'\ 



{n + l + 1/2)1 



(29) 

where L!'n^^'^{r^) is the associated Laguerre polynomial 
and f = r^^. 

The Hamiltonian in Eq. (|28p can be diagonalized using 
a canonical transformation 

N N 
Cnlm^ = E ^L^llrn^ + i'^^ E ^iN+dl^mi (30) 

and 

N N 
i=l J = l 

(31) 

which yields the diagonalized Hamiltonian 



jl'm,a 



(32) 



The index j in ^jima corresponds to the enumeration 
of the quasiparticle states, the index / is the angular mo- 
mentum, and m = —I, — Z -|- 1, . . . , / is the z-component of 
the angular momentum. The index a does not have the 
meaning of a physical (pseudo)spin anymore, but it has 
an auxiliary function. 

The equation for the order parameter profile at zero 
temperature is 



A(r)=5oE 



nn'lj 



21 + 1 
An 



Rniir)Rn'iir)W^^W^ 



j N+n',j^ 



(33) 



which needs to be solved self-consistently together with 
the number equations for the local atom densities 



97 1 1 

(r) - E -r^Rni{r)Rn'iir)WiM',„ (34) 



nn'lj 



and 



nn' Lj 



An 



Rni{r)RMr)WUN.^W,[,+^^^. (35) 
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FIG. 1. (Color online) Order parameter profiles A(r) ob- 
tained using the Bogoliubov-deGennes method for different 
interaction strengths. Open and closed shell data correspond 
to different total atom numbers, A'^ — 4930 and A'^ = 5600, 
respectively. Shown are also the profiles that include the 
Hartree energy shift contribution, calculated for A'^ — 4930. 
Figures show data for kpa = —0.4, kpa — —0.6 and kpa — 
-0.8. 



Fig. [T] shows the order parameter profiles A(r) for differ- 
ent atom numbers, interaction strengths, and also with 
the Hartree energy shift. 



We calculate the coefficients W^ j and the eigenener- 
gies Eji numerically. The Green's function Gy (1, 2) can 
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be expressed via Wi : , a-^ \ ^ f ^nN+i \ r> f \ 

^ "'^ where AJr) = Y.n [ u/i ] Rni{r), 

A^i(r) = E„( Jf^';^ J^n^M, and Pz(cos^i2) = 

2lTTSM=-L^LM(^i'¥'i)^iM(^'2,¥'2) are the Legendre 
polynomials and On is the angle between the vectors ri 
and Yi. Here G is the Green's function in the matrix 



— » ^ 2/ + 1 form G = 1 ^ ^ 

G(ri,r2,f7„) = -V^^P;(cos0i2)x V <^it , 

V The kernel Cikij.L{ri,r3,uj) of Eq. ([26]) is (see Ap- 



(36) pendixH]): 
J 



Afc;i,L(''i,?'3, w) = (2X + 1) ^ 



i Li i2 2Li + 1 2L2 + 1 



, y 4tt 47r 

L1L2 



, x+ ?^f(£^JiLi) -r^F(-£^J2L2) , ^+ "F(--g./iLi) " nF(£;j,L2) 

+ A,,i,,.feA,^^^,^. ^ ^ ^^^^^ ^ ^^^^^ + ^ _ ^^^^^ _ ^^^^^ 



(37) 



Here the occupation numbers are given by the Fermi- 
Dirac equation n^iE) — exp(ff£;)+i at the temperature 



given by 



fcsT = -g. Furthermore, i q g g 



are the Wigner 



3j-symbols. Finally, X%Lu^k = ^J^l^M'^^JiL^A^s) 

and X%L2,ij = A± ;(r3)A±;^^_^.(ri). 

In this section we have derived the kernel of the 
Eq. ([25]) . Cihij,L{ri,r3,Lu), as a function of the coeffi- 
cients W^ j and the eigenenergies Eji, which are to be 
calculated numerically. Thus Eq. can be solved, and 
the resulting density response can be obtained. 



D. Basic definitions for numerical calculations 

In this section we will discuss the parameters used in 
the numerical calculations. As noted in Eq. (fT2)) . the gas 
of atoms of mass m is confined in the harmonic trapping 
potential of a frequency cjt- The system is characterised 
by two units: the unit of energy (the difference between 
neighbouring levels) Hut and the unit of length (the oscil- 
lator length) rose = the numerical calculations 
we use dimensionless values, in the unit system based on 
hujj^ and rose- 

The interaction strength is determined by the two- 
body scattering length a as 1/go — ^^2^ — li^) where 
the position dependent renormalization coefBcient is cal- 
culated in the local density approximation [25| and is 



mkc{r) f K{r) 1 + K{r) 

iOST 



27r2?i2 



1 — K(r) 



(38) 
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where ™: = - \muj?,r\ = m - ^--t' 

K(r) = kp{r)/kc{r), and E^. = ftwx (2iVc + 1) is the cut- 
off energy. We use cut-off energy E^ = 60 huT (for runs 
with the Hartree energy we used Ec = lOOfiwx), which 
proved to be sufficient for studying the density response 
of the rather weakly interacting gas, resulting in at most 
1 % error in the magnitude of the pairing field. The Fermi 
energy is Ep — 2Ahu}T for N — 4930. For the interac- 
tion in the density response in Eq. (j26p we use the above 
renormalized interaction go as in Ref. [20j . 



IV. RESULTS 

In this section we present the numerical results of the 
density response. In Section IIV Al we relate the density 
response with the frequency of collective excitations. In 
Section IIVBI we study the density response function for 
weak interactions, and in Section IIV CI we interpret the 
width of a band of collective excitations as the damping 
rate of the modes. In the Section [TV Dl we show how the 
frequencies of the excitations depend on the interaction 
strength and discuss the interesting effect of merging of 
the pair vibration and collisionless hydrodynamic excita- 
tion branches. In the Section FlV El we discuss the cases 
of closed and open shells (i.e. cases of different numbers 
of atoms) . 
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FIG. 2. (Color online) The density response in the extreme 
weakly interacting limit kpa = —0.22 (A (0) — 0.06 Hljt) 
shows two groups of excitations. One, a prominent one, at 
frequencies tj ~ 2 ut and another (shown enlarged in the in- 
set) at tj « 0.15 ojt. Shown are the full density response A{uj) 
and the single particle response Ao{uj). 



FIG. 3. (Color online) The density response for a slightly 
stronger interaction than in Fig. [5] kpa = —0.50 (A (0) = 
0.78fitiJT), shows the increase in the widths of the collective 
excitation bands (in A{ijj)) and the separation of the single 
particle excitations (in Ao{uj)). 



B. Results for weak interactions 



A. Spectrum of the monopole mode 



Our goal is to study the density response 
-4ij,L(ri,r5,a;), which was defined in Section IIIIBI 
As was discussed in Section IIIBl the peaks in the 
response function yield the frequencies of the collective 
excitations of the system. In order to better understand 
the physical origin of various excitations we will consider 
also the non-RPA response .4oij,L(7'i, rs, w), which was 
discussed in Section IIII Bl The peaks in this function 
reflect the frequencies of single particle excitations. 
Thus in this article we will call it the single particle 
density response, and ^y^L(ri, rs, w) the full density 
response. The simultaneous study of Aij^Liri^r^^uj) 
and AQij^L{ri,r^,uj) allows one to identify the origin 
of different collective excitations; some excitations in 
Aij^L{ri,r5,uj) have their origin in ^oy,L(?'i, "^5, w) as 
corresponding single particle excitations, while others, 
which we refer to as purely collective excitations, do not 
have a corresponding peak in Aoij.L{fi, f5,uj). 

The density response ^ij_i(ri, rs, cj) depends on six 
parameters: spin indices i and j, the angular momentum 
L, positions ri and rs, and the frequency lo. In addition, 
the response can be calculated for three different kinds 
of probing fields (p-f, (pi and r/, as described in Eq. (1151) . 
However, in Section [IV Fl we will show that the main fea- 
tures of the collective excitation spectra do not depend 
on the parameters i, j, ri, nor on the fields 0^, 4>i, rj, 
but only on the angular momentum L and the frequency 
w. In this paper we will study only the monopole mode, 
corresponding to L = 0. Nevertheless, our method allows 
also the calculation of the response function for higher an- 
gular momenta L > 0. Hence we will now limit the study 
to Af^{ri = 0, rs = 0,a;)L ^^g, denoting it as A{uj). 



In this chapter we study the full density response 
^-|-|-(a;) = 1^ and the corresponding single particle den- 
sity response ^o,tt('^)' which we denote as A{u!) and 
Ao{ui), respectively. From here on until further notice 
we consider the response for the total number of atoms 
N = 4930. 

Typical results for the density response in the weakly 
interacting regime are shown in Figs. [2] and [3l Fig. [2] 
shows the full density response A{uj) and the correspond- 
ing single particle density response Ao{uj) for kpa — 
—0.22. This case corresponds to the gap in the centre 
of the trap A (0) = 0.06 Hujt. In Fig. [3] are shown re- 
sponses for a higher interaction fcpa ~ —0.50 (the gap in 
the centre of the trap A (0) = 0.78 Itiut). 

Both figures show two groups of peaks, one at low fre- 
quencies uj < 2ujt and another close to 2 lot. This is a 
typical result for sufficiently weak interactions, with the 
key criterion being the value of the gap A(0) < 2Tiujt. 
The group of peaks in the vicinity of 2 wt is present both 
in the full response A{uj) as well as in the single par- 
ticle response Aq{uj), the corresponding single particle 
excitations describing transitions of single atoms from 
the Fermi surface to the next higher oscillator energy 
level. These excitations are described by the hydrody- 
namic model for a collisionless gas. The other group 
of peaks at lower frequencies is shown in the insets of 
Figs. [2] and [3l These peaks are present only in the full 
density response A{uj) but not in the single particle den- 
sity response Ao{uj) and thus these are purely collective 
excitations. These are called pair vibration modes and 
they describe pair amplitude modulations [l^. In the 
very weakly interacting limit the pair vibration mode fre- 
quency is given by the pairing gap as cj « 2A(0)/fi-, but 
as already seen in Fig. [3] the relation does not hold for 
stronger interactions (see the positions of the A{l^) peaks 
in Figs. [2] and 13] compared to 2A(0)). 

Note that for the pair vibration modes we do not ob- 
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FIG. 4. (Color online) The pair vibration modes, corre- 
sponding to purely collective excitations, for kpa = —0.39 
(A(0) = 0.30 fe^T) and fcpa = -0.67 (A (0) = 2.26 &»jt). 
Shown are also the corresponding collective mode band widths 
0.06 cjt and 0.22 a;T, respectively. 

serve an individual peak but rather a group of peaks close 
to each other. With the interaction increasing, the peaks 
are shifted to higher frequencies and the distance between 
them increases. We discuss the interaction dependence 
more in Section HVDI 

Also the group of peaks at w « 2 experiences sig- 
nificant changes when the interaction is increased. Fig. [3] 
shows how for stronger interactions the three single par- 
ticle excitations are accompanied by at least 17 peaks, 
corresponding to collective modes, between w = 2.01 
and to = 2.78 wt- 



C. Damping rate 

Excitations with similar frequencies can be coupled, re- 
sulting in a damping of the modes. Even though the exci- 
tations for a finite system manifest as poles with real en- 
ergies and vanishing imaginary parts, in practice a large 
number of collective excitations with nearby lying fre- 
quencies cannot be distinguished. The distance between 
the mode frequencies within the band yields the time re- 
quired for resolving the various peaks, but for time scales 
shorter than this the resulting real time evolution ap- 
pears damped. For short time scales, the characteristic 
damping rate is given by the width of the collective ex- 
citation band. Fig. 2] shows the pair vibration modes for 
stronger interactions, revealing the gradual increase in 
the collective excitation band width and thus the increase 
in the corresponding damping rate. For fcpa = —0.39 
(A(0) = CSfix^x) the width is roughly 0.06 and for 
kpa = -0.67 (A(0) = 2.26 /iwt) the width is 0.22 wt- 

Fig. [5] shows the width of the lowest frequency collec- 
tive excitation band (the pair vibration mode for weak 
interactions) as a function of interactions (both as a func- 
tion of fcpa and A(0)). While for weak interactions the 
band is narrow, the width increases rapidly when ap- 
proaching the crossover regime kpa w —0.8, correspond- 
ing to A(0) w AhojT. Close to the crossover, where the 




FIG. 5. (Color online) The lowest frequency collective excita- 
tion band width as a function of interactions. The band width 
increases rapidly when the order parameter A(0) approaches 
4/kjT at kpa ~ —0.8. 

pair vibration and the hydrodynamic modes merge, the 
width of only the lowest excitation band may not be a 
good measure of the damping as the distance between 
the two collective mode bands is less than the widths of 
the two bands. 

The number of peaks within the collective mode band 
depends on the system size, with larger systems yield- 
ing more peaks (the procedure of defining "a peak" is 
described in Section [C|. In the thermodynamic limit 
— )■ CJO, we expect that the collective mode band no 
longer consists of separate peaks but becomes continuous. 
However, as long as the magnitude of the pairing gap is 
unchanged the width of the band is unaffected. This con- 
jecture is supported by our calculations with higher atom 
numbers (up to iV = 49300). 



D. Interaction dependence 

When the interaction strength increases, the two 
branches, the pair vibration modes and the collisionless 
hydrodynamic modes, approach each other and eventu- 
ally merge. Instead of the Fermi energy, the relevant 
energy scale in this merging regime is the trap oscilla- 
tor energy Hlot. Subsequently the interaction strength is 
best measured as A{0)/huJT instead of the standard kpa. 

The main result of this work. Fig. ^ shows how the po- 
sitions of the peaks evolve as a function of the interaction 
(showing both fcpa and A(0)). As discussed above, in 
the weakly interacting regime (here defined as the regime 
where A(0) < 2 Hujt) there are two main branches of ex- 
citations. One branch originates from the collisionless 
hydrodynamic excitation (or single particle excitations) 
at frequency uj = 2 ut and the other branch describes 
pair vibration modes, starting from the zero frequency. 
For weak interactions, the latter follows the frequency 
2A(0)/fi [l^, but the branch deviates from this asymp- 
totic value already at A(0) « 0.2 Hlut. The collisionless 
hydrodynamic and the pair vibration branches merge at 
around fcpa —0.8, corresponding to A(0) ~ ifkuT- For 
interactions beyond that, the hydrodynamic and the pair 
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FIG. 6. (Color online) The peaks in the density responses as a function of the interaction kpa or A(0). The crossover from 
the pair vibration mode to the strongly interacting hydrodynamic mode occurs for A(0) ~ 4&jt when the low frequency pair 
vibration mode merges with the weakly interacting coUisionless hydrodynamic mode. 



vibration modes cannot be distinguished anymore and 
the two modes transform into the strongly interacting 
Goldstone mode in the strongly interacting limit. 

The crossover between the pair vibration mode and 
the Goldstone mode takes place when the order param- 
eter A(0) is of the order of a few trap oscillator energies 
hujT. Since the key energy scale is given by the trap 
frequency instead of the Fermi energy, the crossover is 
realized at different interaction strengths kpa if wt is dif- 
ferent. This was discussed already in Refs. The 



interaction strength studied in Ref. [l^ corresponded to 
A(0) < 2 huT, whereas in Ref. [1^ the monopole mode 
was studied for A(0) « 6 Hujt- Here we consider the 
crossover region between these two limits. 

Notice the depression of the collective mode frequen- 
cies in the crossover regime in Fig. [S] and the increase 
of the collective mode band widths in Fig. [5j These ef- 
fects, the suppression of the collective mode frequency 
and the increase in the damping rate in the crossover 
regime where A(0) ~ Hlot, are in qualitative agreement 
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with similar effects observed in experiments |10l - ll2| al- 
though the symmetries of both the trap and the pertur- 
bations are different in the experimental setups. The ex- 
periment in Ref. lo'] observed a dramatic increase in the 
damping rate and a decrease in the radial compression 
mode frequency at a magnetic field B ^ 910 G corre- 
sponding to the value fcpa « —0.5. The radial breathing 



mode was studied in Ref. |ll| , revealing a strong damping 
rate maximum and a depression in the collective mode 
frequency at a magnetic field B ~ 1080 G, corresponding 
to kpa ~ —0.74. Finally, the radial quadrupolc mode was 
studied in Ref. [l^ , exhibiting a damping rate maximum 
and a depression and a jump of the collective mode fre- 
quency at a magnetic field B w 950 G, corresponding to 
kpa « —0.8. Despite the differences between the exper- 
iments in the critical interaction strengths for observing 
the damping rate maxima, all three experiments are in 
the regime where the predicted BCS pairing gap A{0)/h 
is of the order of a few oscillator frequencies. However, 
a detailed comparison with the experiments would re- 
quire a theoretical study of the same radial modes, not 
accessible in the present model assuming the spherical 
symmetry. 

Fig.[7]shows the effect of the Hartree energy on the den- 
sity response. The lowest energy collective mode band is 
shifted to higher frequencies and it is wider. This re- 
sults in the crossover between the two collective mode 
branches occuring at a weaker interaction than when the 
Hartree energy was neglected. However, the qualitative 
features are unchanged from Fig. [6] 



E. Open and closed shells 

Depending on the number of the atoms, the single par- 
ticle density response Aq can show significantly different 
behaviour. In this section we discuss the cases of closed 
and open shells [l^. The closed shell, which was the 
case discussed above with N = 4930, corresponds to a 
case in which the uppermost occupied energy level at the 
Fermi surface is fully occupied in the ideal noninteracting 
system. In contrast, in an open shell configuration (for 
example N — 5600 used below) only part of the energy 
level at the Fermi surface is occupied. 

The two cases give very different single particle den- 
sity response functions Aq in the weakly interacting limit 
(A(0) < 2 Hujt)- In the case of the closed shell, the lowest 
frequency single particle transition is the transition from 
the Fermi surface to the next higher n-level, correspond- 
ing to the frequency 2 wt • This results in a branch of 
single particle excitations starting at frequency 2 wt in 
Fig. [HI In the case of the open shell, the energy level at 
the Fermi surface is also available for transitions, and the 
lowest energy single particle transition is a pair breaking 
transition in which the n quantum number is not affected. 
The associated energy change is proportional to 2 A(0) in 
the weakly interacting limit, tending to zero for vanishing 
interaction strength. The evolution of the peak positions 



as a function of interactions for the open shell configura- 
tion is shown in the Fig. |S1 revealing one single particle 
excitation branch starting from zero frequency and two 
single particle branches from the frequency 2 wt • While 
the difference between the two shell configurations is a 
mesoscopic effect, the effect on the single particle excita- 
tion spectrum is artificial as the experimentally observ- 
able quantity is not but the full response A. Indeed, 
comparing the full density responses of the two configu- 
rations does not show any apparent difference, albeit the 
data for the open shell configuration is somewhat more 
noisy due to numerical issues. 



F. Other response functions 

In a trapped and nonuniform system the density re- 
sponse depends on the positions at which the perturba- 
tion is applied on and where the response is measured 
at. In this work we have considered the response at the 
centre of the trap, i.e. we defined 



^(w) = ^tt(^i-0,r5 = 0,w)| 



(39) 



M = / dr^drs A^^{ri,r5,u;)\^^j^^^. 



The reason for this choice is the simplicity and the nu- 
merical efficiency, allowing the calculation of plots such 
as Fig. [B] In contrast, the standard measure of the den- 
sity response is the strength function, which is obtained 
by integrating the density response over the whole trap 



(40) 



The two definitions, Eqs. ([39]) and (|40| produce slightly 
different results, but the main features of the full density 
response are the same: the frequencies and the widths 
of the collective mode bands are the same for both as 
shown in Fig. |9l The reason for the similarities between 
two methods is that the generalized random phase ap- 
proximation couples the excitations at the centre of the 
trap to excitations all over the system, and hence one 
can excite for example surface modes even by consider- 
ing only the centre of the trap. Such indirect effects come 
at the cost of reducing the amplitude of the correspond- 
ing collective mode peaks. However, when considering 
only the frequencies and the widths of the collective ex- 
citation bands, the actual amplitudes of various peaks 
are irrelevant. 

In contrast, looking only at single particle excitations 
Aq, the local response (at the centre of the trap) has only 
three low energy branches, whereas the 'single-particle' 
strength function Aq = J drxdr5Ao{ri,rc^,uj) has a large 
number of peaks, corresponding to single-particle excita- 
tions at different parts of the trap. Without the coupling 
provided by the generalized random phase approxima- 
tion, the single-particle excitations localized at the cen- 
tre of the trap will remain localized. The single particle 
excitation spectra are thus clearly different in the two ap- 
proaches. However, the experimentally relevant quantity 
is the full density response. 
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FIG. 7. (Color online) The peaks in the density responses as a function of the interaction fcpa or A(0), including the Hartree 
energy contribution. The Hartree energy compresses the gas, resulting in a larger order parameter A(0). Comparing with 
Fig. [SI this shifts the crossover from the pair vibration mode to the strongly interacting hydrodynamic mode to lower interaction 
strengths. 



The density response results discussed in this 
manuscript have considered only the response to prob- 
ing field directly affecting only the density of spin-f 
atoms. Fig. [TU] shows the response for all three probing 
fields (/()^, (j)]^ and 77, as described in Eq. (|T5|) . While the 
magnitudes of the individual peaks are different for dif- 
ferent probing fields, the main features of the responses 
are the same, namely the frequencies and the widths of 
the collective excitation bands. The same is true also 
when considering the response of the pairing field. 



V. DISCUSSION 

In this work, we have used the random phase approx- 
imation together with the self-consistent Bogoliubov- 
deGennes method for studying the density response of 
a trapped Fermi gas. Concentrating on the monopole 
mode in a spherically symmetric trap we have analysed 
in detail the interesting crossover regime in which the 
coUisionless gas becomes strongly interacting. We ob- 
serve the merging of a pair vibration mode, originating 



from a low frequency excitation in the extreme weakly in- 
teracting gas, and a coUisionless hydrodynamic mode at 
frequency 2a;T, the combined collective mode eventually 
becoming the strongly interacting hydrodynamic Gold- 
stone mode in the unitary regime. The merging of the 
two collective mode branches is signalled by a depression 
of the mode frequencies and an increase in the damp- 
ing rate, in good agreement with the experiments done 
in elongated traps. In the future it will be interesting 
to generalise the method to non-spherically symmetric 
systems [26l - l29j , allowing the study of various radial and 
axial modes and thus a more detailed comparison with 
the reported and possibly new experiments. Other inter- 
esting future extensions of the method are studies of the 
FFLO state [13,^1, sensor applications [s^] and dimen- 
sionality effects |33| . 
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Appendix A: Density response within the Random 
Phase Approximation 

In this appendix we discuss the derivation of the den- 
sity response Eq. (j23l) . starting from the Hamiltonian in 
Eq. dn]). The Green's function G(l,2) satisfies the fol- 
lowing equation 



d 



V? 
2m 



T3 G(l, 2) = 16(1 -2) + VF(1)G(1, 2) + 5oAfint(l, 2), 



(Al) 



where / is the identity operator. 



and 



Wil) = 



-eMl) 77* (1) 
^(1) e(/);(l) 



(A2) 



Mint (1,2) = 

;r^j(2)n(l)^t(l 

;r^j(2)v^|(i)n(i) 



TV';(2)n(l)V'^(l)^ 
TV'I(l)n(l)V;(2) 



(A3) 
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where ( T . 



is a time-ordered correlator. The bare 



Green's function Go(l,2) is defined similarly but in the 
absence of the perturbation and the interactions go = 



d 

dri 



2rn 2 
Now, we can express Eq. (jAip using Go(l, 2) as follows 



T3 ] Go(l,2) =M(l-2). 

(A4) 



G(l,2) = Go(l,2) 

d3 j d4Go(l,3)VF(3),5(3-4)G(4,2) 

.go / d3 / d4Go(l,3)S(3,4)G(4,2), 



(A5) 



where the self-energy E(3, 4) = / d5Afint(3, 5)G-i(5, 4). 
Eq. (|A5p can be formally written as 



G(l,5) ^ Q J G(5,2) for h equal to 0^, 0^ and ry, 
respectively. 

After transforming G(l,3) = T3G(1,3), Eq. ([XT]) be- 
comes: 

Ay(l,2,5) = io^j(l,2,5) 

+ 90 J rf3Lifefej(l, 2, 3)Jl;/(3, 3, 5) j-^g^ 

-5-0 / d3L,fe,j(l,2,3)ife/(3,3,5), 



where i.,(l,2,5) = £,^,,(1,2,3) = 



G,fc(l,3)Gy(3,2) 
10 



G(l,5) 
G(l,5) 





1 

1 



and 

G(5,2), G(l,5) 



Ao(l,2,5) 

1 



G(5,2), 



G(5,2) for h = (j)^ and 77, respec- 



tively. 

Appendix B: Coefficients Cikij.L{ri,rs,,u}) via BdG 
Green's functions 



G-i(l, 2) = G^\l, 2) - W{1)5{1 - 2) - 2). (A6) 



The Fourier transform of the matrix element Likij is 
Likij{vi,Y'i,uj) = ^X^n,, Gjfe(ri,r3,f7„)Gzj(r3,ri,w 

iln), where /3 = is the temperature, and 



The next step is to apply the variational deriva- 



fir. 



— '•^"^^•'^ are Matsubara frequencies. As dis- 
j^gy to both sides of Eq. (jA6|) for separate cussed in Section IIIIBl the spherical symmetry al- 
lows a great simplification of the problem, and one 
needs only to solve the coefficients Cikij,L{i'i,r^,uj) ~ 
/d7i3sin7i3Life;j(ri,r3,iLj)PL(cos7i3), which can be 
written as 



five 



cases h = 



Applying 



-/ d3d4G(l,3)^^^^G(4,2) yields the equation 
<5G(1,2) 



5h{5) 



= ^0(1,2,5) 



^3S(1,3) -.3(-S + ^r3) G(3,2). 



5h{b) 5h{b) 



Cikij,L{ri,r:i,u) = (2L-I- 1) --^ ^ 

L1L2 

Gife,Li(ri,r3, f2„)G;j,L2(r3, ri,tj -I- il„), 



L Li L2 




(A7) 



(Bl) 



where 



Here Ao(l,2,5) depends from the field h, and is 



L Li L2 




are Wigner 3j coefficients. 



G(l,5) 



1 




G(5,2), 



G(l,5) 




1 



G(5,2), 



The Bogoliubov-deGennes Green's functions are pro- 
vided by Eq. (1551) . yielding the final result 



J 



Cikij,L{ri,r3,uj) = {2L + 1) ^ 

L1L2 



L Li L2 




2Li + 1 2L2 + 1 



An 



47r 



J1J2 



+ ^.JlLi,lkKj2L2,lj 



J2L2.13 ^ + ^^^^^^ _ ^^^^^ 



+ '^JiLi,ik^J2L2,lj 



W + Ej^Li + Ej^L2 



E.J2L 



JiLi.ik .72-^2, y U] ~ E 



JiLi 



E 



J2L2 



(B2) 



where A 



Aj^^^ ;(r3)AjJ^^ j(ri), and the quasiparticle wavefunc- 
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tions are defined in the main text below Eq. (l36l) . 



Appendix C: Procedure of defining tlie pealcs 

Numerical calculations with very high precision show 
that the peaks in the density response have the shape 
~ where wo is the (yet unknown) frequency 

of the mode. However, for most of the figures in this 
work, we have limited the frequency grid resolution to 
Suj — 3 ■ IO^^ojt- This resolution is high enough to 
show the existence of a peak but not to show the de- 
tailed shape. For this resolution, some modes result in 
extremely high peaks (if the frequency grid happened to 
coincide with the mode frequency) regardless of the ac- 



tual prefactor or the amplitude. We decided to define 
'a peak' as a frequency for which the derivative of the 
density response is higher than some chosen cut-off. The 
reason for using the derivative instead of the height of 
the peak is that the base level of the response (the value 
of the response in areas away from any peaks) is slowly 
changing with increasing frequency w. Such slow changes 
in the base level do not affect the derivative and the peaks 
are more clearly seen. The caveat is the higher sensitiv- 
ity to numerical noise. For creating Figs. [S] and [5] we 
mark a peak when the derivative > 3 • 10^. The 

choice of this value is a compromise between limiting the 
required resolution and in avoiding the effect of the nu- 
merical noise. The main results, such as the frequencies 
of the collective modes or the widths of the excitation 
bands, are not sensitive to this choice. 



